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ABSTRACT 

This study is concerned with the early evolution of magnetic fields and differential rotation 
of intermediate-mass stars which may evolve into Ap stars. We report on simulations of the 
interplay of differential rotation and magnetic fields, the stability limits and non-linear evolu- 
tion of such configurations, and the prospects of dynamo action from the unstable cases. The 
axisymmetric problem delivers a balance between field amplification and back-reaction of the 
magnetic field on the differential rotation. The non-axisymmetric case involves also the Tayler 
instability of the amplified toroidal fields. We consider limits for field amplification and apply 
these to young A stars. 

Apart from its application to Ap stars, the instability is scrutinized for the fundamental 
possibility of a dynamo. We are not looking for a dynamo as an explanation for the Ap star 
phenomenon. The kinetic helicity is concentrated near the tangent cylinder of the inner sphere 
of the computational domain and is negative in the northern hemisphere. This appears to be 
a ubiquitous effect not special to the Tayler instability. The latter is actually connected with 
a positive current helicity in the bulk of the spherical shell giving rise to a small, but non- 
vanishing a-effect in non-linear evolution of the instability. 
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1 INTRODUCTION 

Stellar objects involve in many cases the amplification of magnetic 
fields by differential rotation. The differential rotation at the bottom 
of the solar convection zone, for example, acts as a generator for 
toroidal magnetic fields while the differential rotation is imposed 
from the convection zone. In principle, there is a large domain of 
radiative zones in stars in which the amplification by differential 
rotation is conceivable, where no convective motions disturb the 
winding-up. This concerns the whole range of intermediate-mass 
and massive stars from 1.2 solar masses up to several tens of a solar 
mass. All these stars possess large radiative envelopes (apart from 
a possible, very thin convective surface layer). We will often speak 
about A and Ap stars in this Paper, but the simulations are actu- 
ally applicable to a wider range of stars possessing non-convective 
domains, such as solar-like stars with radiative interiors or more 
massive B stars. 

According to Stepien (2000), the pre-main-sequence evolu- 
tion of stars with masses near two solar masses is accompanied 
by changes in angular momentum through disk accretion, magnetic 
star-disk coupling, and magnetized winds. All these effects likely 
lead to a differentially rotating interior. When approaching the zero- 
age main sequence (ZAMS), the stars have established their radia- 
tive envelopes, and accretion has ceased. It is then mainly by mag- 
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netized winds that the angular velocity is reduced. This, however, 
affects only the surface, while it takes extremely long for the small 
plasma viscosity to reduce also the angular velocity in the interiors. 
The coupling of the surface to the interior is strongly enhanced by 
magnetic fields. There are a few cases for which an increase of the 
surface rotation period has actually been observed, such as V901 
Ori (Mikulasek et al. 2008). 

Another source for differential rotation, especially for stars 
with masses less than 2.5M0, may be the convective phase in 
the pre-main-sequence evolution. While there is no specific model 
available for a progenitor of A stars, the computations by Kiiker & 
Riidiger (2008) indicate tendencies of stronger differential rotation 
for faster rotation, hotter surface temperatures, and especially thin- 
ner convective envelopes. To which extent the differential rotation 
will be present in the remaining radiative star is debatable, since 
small magnetic fields in the (then) radiative interior are sufficient 
to prevent differential rotation from spreading below the convec- 
tion zone (Riidiger & Kitchatinov 1997; Gough & Mclntyre 1998). 

Studying the interaction between differential rotation and 
magnetic fields is important for the understanding of the existence 
of magnetic Ap stars, and for the explanation why only a fraction 
of A stars show the peculiarities. Differential rotation or magnetic 
fields or combinations of the two may become unstable against 
small perturbations (e.g. Watson 1981; Oilman & Fox 1997; Dik- 
pati et al. 2004; Braithwaite 2006b; Arlt et al. 2007; and many 
others). Where MHD is a suitable description of the astrophysical 
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objects, we can group the instabilities into shear-driven, current- 
driven, and buoyancy-driven instabilities. 

Current-driven instabilities are particularly interesting for 
toroidal magnetic fields, since these typically reach the largest 
strengths if differential rotation is present. Such toroidal fields also 
come along with currents giving rise to a kink-type instability, 
whose most easily excited modes are typically modes with low 
azimuthal wave number (Vandakurov 1972; Tayler 1973; Spruit 
1999). We refer to this instability here as Tayler instability. 

Dynamo action arising from the Tayler instability was sug- 
gested heuristically by Spruit (2002) and later in simulations by 
Braithwaite (2006a) who found sustained, but not constant mag- 
netic field energies in a differentially rotating cylindric domain. 
Spherical, three-dimensional simulations of the solar radiative in- 
terior have been performed by Brun & Zahn (2006) where the evo- 
lution of a fossil field together with an imposed rotation profile at 
the top of the domain (i.e. the bottom of the convection zone) was 
investigated. The initial poloidal magnetic field as well as the later 
combination of toroidal and poloidal fields showed distinct insta- 
bility signatures. The simulations were examined by Zahn, Brun 
& Mathis (2007) on possible dynamo action, but none was found. 
The difference between the two attempts may be the presence of 
an enforced differential rotation inside the domain of Braithwaite's 
(2006a) simulations, as opposed to the enforcement on the bound- 
ary only in the computations by Zahn et al. (2007). Another non- 
linear simulation of the Tayler instability was performed by Gellert, 
Riidiger & Elstner (2008) who also used a cylindrical setup with a 
differential rotation enforced at the outer radial boundary, but in this 
case depending on vertical direction z only. In terms of a dynamo- 
q they find values of \a zz \ ~ 0.05 measured against an imposed 
vertical field. That \a zz | is 100 per cent of the rms velocity imply- 
ing that the flow is entirely helical. However, all magnetic modes 
decay when the external magnetic field is switched off. 

The precise nature of the non-linear evolution of the Tayler 
instability with a possible dynamo effect is thus a very interesting 
issue on which we follow up here. The scenario of amplification of 
magnetic fields by differential rotation, the stability of the result- 
ing configurations, and the possible dynamo effect in case of insta- 
bility is very general to non-convective zones in a wide range of 
stars. Section 2 describes how the computations were set up, while 
Section [3] shows the results of the axisymmetric, non-linear simu- 
lations for the field-amplification. The analysis in Section [4] deals 
with the stability of the results in Section[3] and Section|5]evolves 
the full three-dimensional and non-linear problem. Dynamo coef- 
ficients are derived in Section[6] With respect to the instability, we 
will focus on the pre-main-sequence evolution of Ap stars here, 
whereas the possible dynamo action is not believed to be an im- 
portant Ap star phenomenon, but is of general interest for dynamo 
theory. 



2 NUMERICAL SETUP 

The computational domain is a spherical shell with the radius r run- 
ning from an inner radius n to an outer radius r a which is normal- 
ized to one by the stellar radius R in the following. The colatitude 
8 always covers the range from to n; the azimuthal coordinate <j> 
runs from to 2-7T in the non-axisymmetric simulations. The inner 
radius was set to n = 0.5 in all the computations of this study. 

The investigation involves three steps: Step I one is an axisym- 
metric simulation with an initial differential rotation fits) where 



s = r sin 6 is the axis distance, and an initial poloidal magnetic 
fieldB = (B r (r,9),Be(r,9),0). 

Step II is a stability analysis where a number of snapshots of 
the simulation of the first step is used as background states for a lin- 
ear analysis of the stability of the system against non-axisymmetric 
perturbations. We will restrict ourselves to the stability of perturba- 
tions with azimuthal wave-numbers of m = 1. 

Step III consists of fully non-linear, three-dimensional compu- 
tations of the evolution. In these simulations, the non-axisymmetric 
perturbation is either present from the beginning, or is imposed at a 
later step as an external 'kick' to the system. We can thus compare 
the evolutions of states which are either Tayler-stable or Tayler- 
unstable according to Step II. The non-linear evolutions also al- 
low to assess the possibility of driving a dynamo with the insta- 
bility. Dynamo effect will be visible only in terms of non-vanishing 
dynamo-coefficients of a mean-field description, but not in growing 
and sustained magnetic fields, because our system lacks a source of 
energy. 

The computations employ the spherical MHD code by Holler- 
bach (2000) which integrates the momentum, induction, and tem- 
perature equations in the Boussinesq approximation. We are not 
using temperature fluctuations in this study; the normalized, incom- 
pressible MHD equations are thus 

-(« ■ V)u + (V x B) x B 
PmAii, 



dt 



-Vp 



^ = V x (« x B) 



AS, 



(1) 
(2) 



where u and B are the velocity and magnetic fields and p is the 
pressure. Additionally, the relations V ■ u = and V ■ b = 
hold. The magnetic permeability and the density are set to unity in 
this system of units. Since all quantities are normalized with the 
magnetic diffusivity r\, the magnetic Prandtl number Pm = u/rj 
appears in the momentum equation, where v is the kinematic vis- 
cosity. Velocities are normalized by rj/R, times by R 2 /r/. 

For Steps I and III, the differential rotation and the poloidal 
magnetic field are initial conditions. They evolve freely without 
any further imposed properties. The profile of the initial angular 
velocity follows 

where Q,q is the non-dimensional angular velocity on the axis and q 
is a parameter controlling the steepness of the rotation profile. With 
our normalization, the initial magnetic Reynolds number is actually 
Rm = Oq, since 



Rm = 



Tj 



(4) 



where SI* is the angular velocity of the star in physical units. The 
magnetic Reynolds number was Rm = 20 000 in all computations 
except one; the magnetic Prandtl number was always Pm = 1. 

The pre-main-sequence evolution of the rotation period of 
intermediate-mass stars contains most likely both a spin-up (con- 
traction and accretion) phase and a spin-down phase with magne- 
tized winds (Stepieri 2000). We are concerned with the evolution 
near the end of the pre-main-sequence phase and the early main- 
sequence life when most of the star is already radiative (convec- 
tively stable). The probable angular-momentum change is a spin- 
down then, acting on the surface of the star, and this is the reason 
why we assumed an angular velocity decreasing with axis distance. 

The initial profile uses q(t = 0) = 4, but q will later be 
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used to measure the actual steepness of the rotation curve in the 
simulation. The specific angular momentum s 2 Q(s) is increasing 
with axis distance everywhere in the computational domain. The 
system is thus not prone to the hydrodynamic Rayleigh instabil- 
ity which we do not want to investigate here. The hydrodynamic 
stability has been tested numerically with non-axisymmetric per- 
turbations which were symmetric and anti-symmetric with respect 
to the equator. All situations delivered a decay of the perturbations. 

Magnetic fields are measured in terms of Lundquist numbers, 
which is the same as the non-dimensional Alfven velocity in our 
system of units, 



B = 



RB, 



phys 



(5) 



The magnetic diffusivity 77 in time-dependent simulations typi- 
cally represents a value between the microscopic diffusivity of the 
plasma and the turbulent diffusivity resulting from, e.g., averaged 
convective motions. Of the quantities entering @, 77 is the one 
which is known least. It is therefore best to eliminate r\ by Rm and 
thus retrieve the physical magnetic fields by comparing its Alfven 
speed with the rotational velocity, 



(6) 



Considering a 2M star with a radius of R t — 1.5Rq, a rotation 
period of 10 days (fiphys = 7.3-10 -6 ) and a density of 0.015 g/cm 3 
at 0.75-R, , a magnetic field of B = 100 converts to 1660 G in 
physical units with the magnetic Reynolds number of 20 000 used 
in most of our simulations. 

The boundary conditions for all runs presented here employ 
stress-free conditions for the flow and vacuum conditions for the 
magnetic field at both the inner and outer radii, n and r a . Although 
the inner sphere is supposed to be highly conducting, a vacuum 
condition helps in obtaining smoother solutions compared to a su- 
perconducting condition. The latter prohibits a penetration of the 
magnetic field through the inner boundary. This typically leads to 
strong currents near the inner boundary which are not there if the 
conductivity would change continuously into the inner sphere. If, 
for efficiency reasons, the whole star cannot be computed, the vac- 
uum condition is a suitable choice despite the non-vanishing con- 
ductivity. We will come back to a test run with a perfect-conductor 
boundary at r; in Section[5] 

Note that the assumption of an initially independent angular 
velocity of the vertical axis z — r cos 6 (Tayler-Proudman state) is 
not compatible with a stress-free condition on a spherical surface. 
This causes small meridional flows and results in a small deviation 
from the Tayler-Proudman state. However, this is a much better 
choice than an initial rotation profile depending on r which causes a 
severe redistribution of angular momentum towards a configuration 
which is very near a Tayler-Proudman state during the first rotations 
of the system. As we did not want to obtain a mixed evolution of 
this hydrodynamic phenomenon with the magnetic phenomena, we 
chose fi(s) instead of more complicated profiles Q.(r, 9). 

The spin-down by winds during the late pre-main-sequence 
phase of intermediate-mass stars is probably not causing a pre- 
cise Q(s) profile, but for the sake of physical clarity, we choose 
a Tayler-Proudman state as an initial condition for the runs. Note 
that the removal of angular momentum by a magnetized wind is not 
implemented in our setup, so the simulations can be understood as 
mimicking a period around the time when the star enters the main 
sequence. 
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Figure 1. Amplification of the toroidal magnetic field under the influence 
of differential rotation, starting from a poloidal field. The solid line refers 
to an initial maximum magnetic field strength of 300, the dotted line is the 
evolution of an initial field strength of 150. 



3 AXISYMMETRIC EVOLUTION (STEP I) 

First we study the amplification of magnetic fields by shear as well 
as the back-reaction of the fields on the differential rotation. The 
two-dimensional, non-linear simulations thus start with an initial 
differential rotation and a purely poloidal magnetic field in a radia- 
tive stellar zone. These initial conditions evolve freely without any 
imposed flows or fields, neither in the bulk of the computational 
domain nor at the boundaries. 

The early phase of the simulation shows a generation and steep 
amplification of toroidal magnetic field through differential rota- 
tion. The generated B r B^, and, to a smaller extent, also the BqB^ 
impose a Lorentz force to the rotational velocity and redistributes 
angular momentum. This is why at the same time of field amplifi- 
cation, the differential rotation starts to decrease, and the toroidal- 
field growth is thus limited. The whole process reaches a maximum 
of magnetic field energy after a few rotational periods. 

Fig. Q] shows the maximum magnetic field in the spherical 
shell as a function of time. The magnetic Reynolds number was 
Rm = 20 000 and the magnetic Prandtl number was Pm = 1. 
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Two runs with initial, purely poloidal magnetic fields with maxi- 
mum strengths of Bo = 300 and Bo = 150 are shown. The lower 
panel of the Figure shows the decay of the differential rotation. We 
assume that the differential rotation roughly follows a profile {3} 
throughout the simulation, with varying parameters fio and q. The 
values of fl(s) in the equatorial plane of each snapshot of the sim- 
ulation are used to fit a best profile <[3j delivering two time series 
for flo(i) and q(t). It is the temporal evolution of q which is shown 
in the lower panel of Fig. Q] Interestingly, the run with a stronger 
initial poloidal magnetic field shows a reversal of the differential 
rotation profile after 0.006 diffusion times. Only after a very long 
period, not plotted here, the shear converges to zero in the entire 
computational domain. 

Maximum internal field strengths are about 24 kG in the run 
with Bo = 300 and about 16 kG in the run with Bo = 100, using 
the conversion {6]l. Much higher magnetic Reynolds numbers en- 
countered in real stars lead to stronger toroidal fields. A run with 
Rm = 40 000 later described in Section[5]reached about 46 kG. 

The amplification time results from the comparison of the 
Lorentz force acting on the differential rotation with the genera- 
tion of toroidal magnetic field by the shear. The linearized action of 
the Lorentz force may be approximated by Af2 ~ AtB r B^ / ppR 
and the linear winding up of toroidal field in the induction equa- 
tion can be estimated by B^ ~ AtAQ.B r /R. Plugging the 
into the first relation, delivers At = R^fpp/ B r . The amplification 
time is independent of the rotation rate, the amplitude of the dif- 
ferential rotation, and the diffusivity. In our dimensionless system, 
the amplification time is thus At = 0.0033 for the initial field of 
Bo = 300. This is very close the the peak time of the solid line in 
Fig. [T] The dashed curve from half the initial field strength peaks 
at roughly twice that period - a bit earlier since diffusion and hy- 
drodynamic redistribution of angular momentum are also at play in 
this numerical setup. There is a slight tendency of smaller magnetic 
Reynolds numbers producing shorter times to reach the peak mag- 
netic field. This is solely due to the non- vanishing viscosity, being 
present in our system of finite Rm, which adds to the reduction of 
the differential rotation and limits the growth eventually. 

A real system with stellar microscopic diffusivities of say 
r\ = 1000 cm 2 /s will show a nearly stationary phase after the am- 
plification, since the Ohmic decay times are of order gigayears. A 
time-dependent numerical simulation cannot model the very low 
diffusivity at the true rotation rate, and the solution is subject to 
a much larger Ohmic dissipation after the maximum toroidal field 
is reached. The following process is thus characterized by a decay 
of both magnetic fields and differential rotation. Since the simula- 
tions preserve angular momentum parallel to the rotation axis, the 
asymptotic state is field-free and with uniform rotation. We are not 
interested in this asymptotic behaviour. Intermediate-mass stars are 
not living long enough to show this final state during their evolution 
on the main sequence. 



4 NON-AXISYMMETRIC, LINEAR STABILITY (STEP II) 

Differential rotation and rotation in general stabilize any current- 
driven instabilities (Tayler instability) as was shown by Riidiger 
& Kitchatinov (2010). Differential rotation together with diffusion 
has a destructive effect on any non-axisymmetric structures. At 
some point during the amplification of magnetic fields by shear, the 
axisymmetric solution of Section [3] may become unstable against 
non-axisymmetric perturbations, as soon as the stabilizing effect of 
the differential rotation ceases. 
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Figure 2. Stability limits for the magnetic fields and velocity fields of each 
snapshot of the axisymmetric evolution of the Bq = 300 case. The solid 
line is copied from Fig. [TJ The stability limit is plotted for a perturbation 
of m = 1, the dotted line refers to a perturbation whose velocity field 
is symmetric with respect to the equator, the dash-dotted line shows the 
antisymmetric perturbation. 



Each snapshot of the run with Bo = 300 from Section [3] was 
thus tested individually on whether or not a non-axisymmetric per- 
turbation can grow on the particular axisymmetric configuration. 
However, from the snapshots, only the azimuthal velocity and the 
toroidal magnetic field, u$ and B^, are used for the stability analy- 
sis. The full problem was found to be numerically demanding and 
is postponed to a future study. The final conclusions are not affected 
by this simplification. 

We consider the linearized equations for a given azimuthal 
wave number m. The snapshots delivered background states in u 
and B which affect the perturbation (u' , B') and lead to exponen- 
tial growth or decay of the perturbation. 

The stability analysis is based on the linear, normalized MHD 
equations 



du' 



dB' 



[u'xVxu + uxVxi/ — V(u' ■ it) 

+S[(V x B') xB + (Vxi))xB'] 
-Vp + PmAu', 

V x (u x B' + Su' x B) - AB' , 



(7) 
(8) 



again with V • u' = and V • B' = 0. The constant density p and 
the magnetic permeability p are also unity in this normalization. 5* 
is a factor which scales the total background magnetic field taken 
from the axisymmetric run. The critical S for instability is deter- 
mined. The background is unstable, if S > 1. Note that we do not 
evolve the system actually in time: the axisymmetric background 
velocity and magnetic fields are constant. The time-dependence ap- 
pears here solely for the test on growing or decaying perturbations. 

The dotted and dot-dashed lines in Fig. [2] are the results of 
the stability analysis. The evolution derived from the axisymmetric 
simulation starting with an initial maximum field strength of 300 
is taken from Fig. Q] and is plotted as a solid line. By S 1 we re- 
fer to a velocity perturbation u' which is symmetric with respect 
to the equator and has m = 1, Al is the corresponding antisym- 
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Table 1. Four non-linear simulations with their parameters as discussed in 
Sections[5]and[6] The truncations for the Chebyshev, Legendre, and Fourier 
modes are given as K, L, and M, respectively. The column t pcrt gives the 
time when the non-axisymmetric perturbation was injected into the system, 
expressed in diffusion times. 



Run 


Rm 


Pm 


K x L x M 


^pcrt 


NL000 


20000 


1 


35 x 80 x 80 





NL003 


20000 


1 


40 x 60 x 60 


0.003 


NL003h 


40000 
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Figure 3. Evolution of a non-axisymmetric perturbation added to the axi- 
symmetric run of Fig.[T]at t = 0. 

metric perturbation. If at any given time both stability lines are 
above the solid line, the corresponding snapshot is stable against 
m = 1 perturbations. The first snapshots with relatively weak 
toroidal fields are all stable. The stability lines cross the solid one 
at about t — 0.0023. There is a minimum marginal stability at 
t = 0.0055 for the SI mode, and at t = 0.0070 for the Al mode. 

After these times, the stability limits increase again. The sta- 
bilization may be due to the change of sign in the shear. The differ- 
ential rotation is then non-vanishing again and may impose a sta- 
bilizing effect on perturbations despite its relatively small, positive 
amplitude. 

Higher m require similar magnetic fields for instability as was 
shown in test runs with m — 2. 
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Figure 4. Evolution of a non-axisymmetric perturbation added to the axi- 
symmetric run of Fig.[T]at t = 0.003. 
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Figure 5. Evolution of a non-axisymmetric perturbation added to the axi- 
symmetric run of Fig.[T]at t = 0.005. 



5 NON-LINEAR EVOLUTION IN 3D (STEP III) 

The initial conditions are the same as for Section [3] but the sys- 
tem is now extended to include also a Fourier decomposition of 
u and B in azimuthal direction. The main difference is a non- 
axisymmetric perturbation hitting the system at a given time to. 
The perturbation is applied to B, has an azimuthal wave number of 
m = 1 and is symmetric with respect to the equator, i.e. B' r (8) = 
B' r {n - 6), B' e {6) = —B'gfjr — 6), and B'^9) = B'^ty - 9). The 
resulting flow is thus antisymmetric and can be compared with the 
A 1 -mode in the linear stability diagram in Fig. [2] A small overview 
of the models is given in Tabled 

Figs. [5}{5] show the evolution of the energies of the lowest 
azimuthal modes, m — and m — 1, for three different in- 
stances of perturbations, t = 0, t = 0.003, and t = 0.005. The 
main feature of the run with the perturbation at the beginning is 
the quicker decay of the energy in the non-axisymmetric magnetic 
field as compared with the other two simulations. Their perturba- 
tions at t — 0.003, and t — 0.005 are within the unstable window 
(cf. Fig. O and show a much more persistent presence in the sys- 
tem. The energy of the magnetic m = 1 mode may even become 
nearly as large as the magnetic m — energy. It is interesting to 
note that the energy in the m = 1 mode of the velocity field decays 
more rapidly than the magnetic energy in both unstable cases. 

Typical horizontal slices of the velocity and magnetic fields 
are shown in Fig. [6] They are plotted for r = 0.75 which is 
half-way between the inner and outer radial boundaries. Azimuthal 
mode numbers m > 1 are obviously at play, but we cannot talk 
about a turbulent state here. The m = 2 contains roughly half the 
magnetic energy of the m = 1 mode throughout the simulation. 
An example spectrum for the kinetic and magnetic energies of all 
individual m-modes is shown in Fig. [7] The energy in the m = 4 
mode is already two orders of magnitude smaller than the energy in 
m = 1. The total contrast between m = 1 and m = 60 is 10 10 in 
the kinetic energy and one order of magnitude higher in the mag- 
netic energy. Strongest magnetic fields appear to be concentrated 
in low latitudes. Since the axisymmetric parts of u and B were 
subtracted before plotting Fig. [6] the surface maps also reflect the 
energy ratio of magnetic m = 1 energy to kinetic m — 1 energy. 
While the ratio of \B r \ to |w r | is 5.3, the energy ratio would be al- 
most 30, being compatible with the m = 1 energy ratio in Fig. [4] 
between t = 0.003 and t = 0.005. The magnetic fluctuations are 
dominating over the velocity fluctuations in these simulations. 

The impact of the non-axisymmetric instability on the rota- 
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Figure 7. Distribution of kinetic energy (solid line) and magnetic energy 
(dotted line) over the azimuthal wave number m for the run NL003 at t = 
0.005, i.e. 0.002 diffusion times after the perturbation was injected. 
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Figure 8. Variation of the steepness of the rotation profile with time: the 
solid line is the result from the axisymmetric run and is the same curve as 
in Fig.[T| the dashed line is the result from the non-linear run with a non- 
axisymmetric perturbation at t = 0.003 = 9.5P ro t, and the dotted line is 
for a perturbation at t = 0.005 = 15.9P ro t. 

tion profile is shown in Fig. [8] An enhanced angular momentum 
transport to reduce the differential rotation is only notable for the 
simulation with a perturbation at t = 0.003. The most obvious dif- 
ference between the unperturbed run and the unstable runs is the 
time it takes to reach fully uniform rotation. Both perturbed simu- 
lations reach a state of uniform rotation after 0.017 diffusion times, 
while it takes the axisymmetric run about 0.05 diffusion times to be 
similarly uniform in rotation (outside the plotting window). 



6 DYNAMO EFFECT 

The following section is not directly related to the pre-main- 
sequence evolution of Ap stars. We are interested in a possible dy- 
namo effect arising from the instability, since this will open very 
interesting possibilities of magnetic-field generation in various con- 
texts. We are not proposing the dynamo effect to be overly impor- 
tant for Ap stars. 

The non-linear interaction of the unstable m = 1 mode may 
generate an axisymmetric amplification effect of poloidal fields out 
of toroidal fields. In that case, the Tayler instability would provide a 



closed dynamo loop besides the generation of toroidal fields by dif- 
ferential rotation. Sustained dynamo action is conceivable as long 
as there is a mechanism that sustains differential rotation, whence 
toroidal fields prone to instability. Rotating convection zones are 
such providers of differential rotation, for example. 

It is assumed here that the dynamo action can be described 
by a non-vanishing turbulent electromotive force (EMF) deliv- 
ering a large-scale axisymmetric magnetic field from the non- 
axisymmetric instability state (which is not necessarily a turbu- 
lent state). Everything axisymmetric is considered large-scale, all 
non-axisymmetric contributions are considered small-scale in this 
context. We are also assuming that the mean-field coefficients only 
relate the large-scale magnetic field B and its first derivatives with 
the EMF, and use the decomposition 

EMF = a~B + 7 x B - /3(V x B) 

-ix(VxB)- K (VB) (sym) (9) 

into symmetric a- and /3-tensors, the vectors 7 and d, and the 
third-rank tensor k acting on the symmetric part of the tensor gra- 
dient of B. 

When restricting to the diagonal elements of a, one notes 
that only a,^ can cause an axisymmetric growth of the large-scale 
poloidal magnetic field from an axisymmetric large-scale toroidal 
field. Remember that the non-axisymmetric Tayler instability is 
now supposed to be hidden in the mean-field coefficients, and is 
not explicitely visible in this description. As long as we are look- 
ing for the driving of an axisymmetric dynamo, the existence of 
the Tayler instability and its properties are supposed to be entirely 
comprised by the mean-field coefficients. 

6.1 Existence of an 

Our first attempt to look for dynamo action thus consists of search- 
ing for a correlation between the 0-component of the turbulent elec- 
tromotive force, tEMF = (u' x B')$, and the ^-component of 
the large-scale magnetic field, (Bj,). The brackets refer to suitable 
averages in space and, ideally, also in time. We simplify the analy- 
sis further by using the azimuthal direction for the spatial average 
and discard time averages. This type of averaging will be denoted 
by overbars, hence tEMF^ = (u' x B'), and B^. The angular 
brackets will be used for averages over entire hemispheres of the 
computational domain in the following. 

Instead of averaging over time, we rather look at the temporal 
evolution of the correlation between tEMF^ and B^. The fluctuat- 
ing quantities are thus derived by u' = u — u and B' = B — B, 
where overbars are always ^-averages. The values of the average 
turbulent EMF and B^, deliver scatter plots from all the radial and 
latitudinal grid points of an entire hemisphere. The slope of the 
regression line as to represent and the correlation coefficient 
were derived from these sets of pairs {tEMF^(r, 6), B#(r, 0)} 
from each snapshot in time and each individual simulation. 

We also computed the kinetic helicity in the two hemispheres 
(N and S) from the velocity fluctuations by 

Hn = (u ■ curfu'^> N , Hs = (u ■ curlu'^ g , (10) 

where averages are taken over entire hemispheres. It will be inter- 
esting to see to which extent the kinetic helicity is related to what 
we measure as an a-effect. Turbulence driven by convection leads 
to an 

a — -\t C oi{u' ■ curki'), (11) 
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Figure 6. non-axisymmetric velocity field (top) and magnetic field (bottom) at r = 0.75 and t = 0.005 from run NL003. The axisymmetric pails (fast 
rotation and strong toroidal magnetic field) were subtracted from the snapshot before plotting. The contours represent the radial components, while the arrows 
are the projected field directions on the (8, <^)-surface. Note that the actual resolution for the computation of the non-linear terms in both horizontal directions, 
9 and </>, is twice as high as plotted here. 



where r cor is the correlation time. This has also been shown in var- 
ious simulations (e.g. Giesecke, Ziegler & Rtidiger 2005; Kapyla, 
Korpi & Brandenburg 2009). Finally, the root mean square (rms) 
values of the velocity fluctuations are computed for the entire com- 
putational domain for the same snapshots in time. 

We will again refer to the three simulations with perturbations 
at t = 0, t = 0.003, and t = 0.005. We recall that in the first 
one, the perturbation was added to a linearly stable state while in 
the other simulations, the perturbations were added to linearly un- 
stable cases. In the last case, the toroidal field was already past its 
maximum value of — 1434, but had still a supercritical strength 
of B$ = 1361 (we recall that the maximum initial poloidal field 
strength was Bo = 300 for comparison). 

Fig. [5] shows the estimate of for the run NL000 with a 
non-axisymmetric perturbation at t — when no toroidal fields 
were present initially. The evolution of this estimate of is com- 



pared with measurements of the kinetic helicity and the rms-values 
of the velocity components. The initial ot^ is solely due to the ini- 
tial perturbation interacting with the system after the first time-step. 
The 0.^ dies out very quickly - in about one rotation period - and 
oscillates around zero with negligible amplitude. The values of 
are then about 200 times smaller than the rms velocity. The correla- 
tion coefficient of the two quantities tEMF^ and B^ in the second 
panel of Fig.[9]is very close to zero all the time. Open symbols are 
actually those cases for which the hypothesis of uncorrelated quan- 
tities, i.e. no Q-effect, holds statistically. The correlation may be 
still significant for the filled triangles. The limit was set such that 
there is a remaining probability of 1 per cent for the filled symbols 
to represent uncorrelated quantities. 

The helicity measured compares well with the rms velocity of 
about 50 and the assumption that the length-scale for the vorticity is 
roughly r Q — n, whence {u 12 }/ (r — 7i) = 5000. The length-scale 
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Figure 9. Dynamo-effect (top-panel) expressed as a^j, measured from the 
simulation with a perturbation at t = 0. Only an initial period of the full run 
of 0.015 diffusion times length is shown here to make the very short pres- 
ence of non-zero visible. The second panel shows the correlation co- 
efficient between the turbulent electromotive force and B^, the third panel 
shows the kinetic helicity and the bottom panel the rms values of the veloc- 
ity fluctuations. 



Figure 10. Dynamo-effect (top-panel) expressed as a^j, measured from the 
simulation with a perturbation at t = 0.003 = 9.5P ro t- Again, the second 
panel is the correlation coefficient, the middle panel shows the kinetic he- 
licity and the bottom panel the rms values of the velocity fluctuations. Note 
that the period covered by these plots is longer than in Fig. [9] and covers 
about 67 rotation periods. 
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Figure 11. Dynamo-effect (top-panel) expressed as a^^ measured from 
the simulation with a perturbation at t = 0.005 = 15.9P r ot- Again, the 
middle panel shows the kinetic helicity and the bottom panel the rms values 
of the velocity fluctuations. 



is probably slightly smaller than the thickness of the spherical shell. 
We also note that the non-axisymmetric motions are dominated by 
their horizontal components. 

The second set of plots in Fig.[l0]is from model NL003 with 
a perturbation at t = 0.003 which is a bit more than two rota- 
tion periods after the magnetic field has reached its Tayler-unstable 
strength. The values of are about 10 times stronger than in 
model NL000, and they are actually increasing during a short pe- 
riod of about one rotation. The helicity is more than 10 times 
stronger than in model NL000 and has the same sign as un- 
til t ~ 0.0075. The correlation coefficient reaches values of —0.6. 
There is apparently a second phase of a small, but significant 
which has the opposite sign than the kinetic helicity. 

The third set of plots in Fig. Qj] shows the simulation NL005 
with a perturbation at t — 0.005 which is at a time when the axi- 
symmetric toroidal field has a strength of B^, — 1361 and is al- 
ready decreasing, and the average differential rotation is already 
near zero. At that time, the linear stability limit is lowest in the 
whole period considered here. The values of are a bit lower 
than those of NL003, but are also first increasing during about one 
rotation period. Again, the correlation coefficient reaches values 
around —0.6 in this run. The sign change, which was also present 
in the of NL003, takes place at t — 0.0088. This is relatively 
earlier than in NL003, when measured from the perturbation time. 

The radial distribution of the a-effect can tell us something 
about the origin of the dynamo action. If one uses only the values 
of {u x B'}^ and B^ for a given r from all the ^-locations in 
a hemisphere, the resulting correlation will be a local one in ra- 
dius. We should bear in mind though, that the become statisti- 
cally less significant, since we obtain a regression line from 47 grid 
points only, in the cases of NL003 and NL005. We plot the corre- 
sponding distributions of estimates versus radius for NL003 in 
Fig-EHwifh a perturbation at t = 0.003. 

Non-zero are mostly found in the inner half of the spher- 
ical shell, at r < 0.75. Because of the vacuum condition on the 
radial boundaries, we do not expect the turbulent EMF to be zero, 
and the derived from the regression method thus need not van- 
ish. A run with perfect conductor boundary conditions at r; was 
performed to evaluate the influence of the inner boundary. The 
is then indeed zero at the inner boundary. The maximum values of 
a^tp near r — 0.6 are about ±90 as compared to about ±110 in 
NL003. We will show spatial distributions of a-components again 
below, obtained with another method. 

Note that there is no point in looking at the total magnetic en- 
ergies in these simulations. Since neither the differential rotation 
nor the initial magnetic field are imposed anywhere in the compu- 
tational domain, there is no long-lived energy source in the system 
which could drive a dynamo noticeable in sustained magnetic ener- 
gies. The concern of this Paper is rather the evolution of stars into 
objects with apparently stationary magnetic fields, such as mag- 
netic intermediate-mass stars. We are not focusing on sustained dy- 
namo action here. 

The spatial distribution of the kinetic helicity and the current 
helicity is shown in Fig. [13] The concentration of the helicity near 
the inner boundary resembles the result from the a-measurement 
which also showed larger values in the inner part of the computa- 
tional domain. More precisely, we can see that the helicity is actu- 
ally concentrated near the place where the tangent cylinder touches 
the inner sphere. The right panel is derived from the run NL003h 
which has the same parameters as the run NL003, however, with 
Rm = 40 000. While the picture is more noisy than the one for 
lower Rm, we also see the the helicity concentration becomes thin- 
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Figure 13. Left: kinetic helicity in the (r, 0)-plane average over the azimuthal direction taken from the run NL003. Middle: current helicity obtained in the 
same way for NL003. Right: kinetic helicity of the run NL003h with doubled magnetic Reynolds number. 



ner with respect to the axis distance. Note that the perturbation of 
this run actually kicks in at a different stage of super-critically as 
compared to NL003, because rotation alters the stability limits. In 
general, however, we have to conclude that a considerable part of 
the kinetic helicity present in the system is due to the geometrical 
setup, namely the presence of an inner cylinder. The situation is 
unchanged in the run with perfect conductor boundary conditions 
at n. 

This is different for the current helicity. A considerable 
amount of positive current helicity is measured in the bulk of the 
northern hemisphere. We believe that it is mostly this positive cur- 
rent helicity contributing to the small, but positive a^-effect in 
the bulk of the northern hemisphere. The regression method is a bit 
crude at this point; the below results from the test-field method will 
show this more clearly. 

It is interesting to note that the sign of and the kinetic 
helicity is the same in the inner half of the radial extent and es- 
pecially near the tangent cylinder. This contradiction with (TlJ I is 
apparently caused by a one-cell helical motion, rather than by an 
average of turbulent motions. Interestingly, Reshetnyak (2006) also 
finds negative kinetic helicity along the tangent cylinder, albeit it is 
a convective, non-magnetic simulation. The feature is more obvi- 
ous in the run with low Rayleigh number, in which convection is 
not reaching a turbulent state. We think that the negative helicity 
near the tangent cylinder of the northern hemisphere is an inner- 
boundary effect, neither related to convective nor Tayler instability 
turbulence. 



6.2 Mean-field coefficients from the test-field method 

The second attempt to look for dynamo action consists of the test- 
field method as described by Schrinner et al. (2007). In parallel to 



the actual non-linear simulation, several equations for a magnetic- 
field like quantity similar to the induction equation are evolved, 

following the action of the large-scale and small-scale velocity, u 

^ 

and u , from the real simulation on given test-fields B T . Because 
of the similarity of the computational setup use here and of the 
simulations used by Schrinner et al. (2007), we implemented the 
test-field method precisely in the same way. Note the different def- 
inition of the signs of a and -y in ([9jl though. 

The separation of large-scale and small-scale contributions to 
u and is again done by azimuthal averages, leading to a large-scale 
field depending on r and 8 only. Nine different test fields give rise to 
27 test-field equations, delivering nine mean electromotive forces 
£ w from u' and the b (l \ Replacing EMF and B by £ {i) and 

Bx' m {9} leads to a system of equations which can be solved for 
the mean-field coefficients. 

The values of the symmetric a-tensor averaged over time are 
shown in Fig. [14] taken from the run NL003 and the period from 
t = 0.0035 to t = 0.0040. The initial phase of the onset of the in- 
stability has thus been omitted. We see that it is a rr which has the 
strongest peak values, followed by agg and a^. In an axisymmet- 
ric mean-field dynamo, a rr and agg are generators of the toroidal 
magnetic field. While a.^ has the same sign as the kinetic helic- 
ity near the inner boundary, the other two diagonal elements have 
the opposite sign. The anisotropy of the a-effect is due to the rel- 
atively strong rotation of the system. The numbers can be directly 
compared to the velocity measurements shown in Fig. II II 

The contributions from the derivatives of B are rather small. 
The coefficients f3 rr , f3gg, and are all positive and do not ex- 
ceed the order or unity. Interestingly, the contributions from several 
pairs of coefficients nearly cancel in the construction of the turbu- 
lent EMF. In particular, these are f3 r( p cancelling with —Sg, and f3g$ 
with S r . 



Amplification and stability of magnetic fields and dynamo effect in young A stars 1 1 



alpha_rr 



alpha_tt 



alpha_pp 




0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 



1.0 



0.5 



0.0 



-0.5 



-1.0 



alpha_rt 



alpha_rp 



alpha_tp 



i 



0.0 0.2 0.4 0.6 0.8 1.0 




0.0 0.2 0.4 0.6 0.8 1.0 



0.0 0.2 0.4 0.6 0.8 1.0 



Figure 14. Coefficients of the symmetric a-tensor obtained by the test-field method. Units are in the dimensionless units of velocity of the simulation. The 
cross-sections shown here are averaged within t = 0.0035 and t = 0.0040 diffusion times. 



6.3 Comparison 

Both methods have their disadvantages. On the one hand, the re- 
gression method for the determination of ignores the possible 
influence of other mean-field coefficients on the 0-component of 
the turbulent EMF. The test-field method, on the other hand, relies 
on turbulence which is not magnetically driven. Only the veloc- 
ity field from the actual simulation enters the determination of the 
mean-field coefficients. Our simulations are concerned with mag- 
netically generated fluctuations u' and B', however. It is therefore 
very likely that the mean-field coefficients determined by the test- 
field method are systematically underestimated. A comparison of 
the rms magnetic fields, B rms = \J (B' 2 ), with rms velocities, 
?irms = y (it' 2 }, over time reveals a factor of about three between 
the two for the run NL003. 

It is elucidating to compare the turbulent EMF measured from 
the simulation directly with the turbulent EMF constructed from the 
mean-field coefficients together with the measured large-scale B. 
The directly measured EMF is about four times larger than the re- 
constructed EMF. This ratio is very similar to the above mentioned 
ratio of B vmB to w rms . The current helicity shown in Fig.[T3]in the 



shell from r = 0.75 to r = 0.95 is even ten times larger than the 
kinetic helicity in the same region. 



7 SUMMARY AND DISCUSSION 

We studied the amplification and stability of magnetic fields in 
the pre-main-sequence evolution of intermediate-mass stars. Espe- 
cially the phase when the stars possess an extensive radiative en- 
velope already is considered. The amplification of a poloidal mag- 
netic field by differential rotation is studied in a non-convective, 
spherical shell. The feedback by the Lorentz force diminishes the 
differential rotation brings the field growth to a halt. The amplifica- 
tion time and amplitude only depend on the initial, poloidal mag- 
netic field. The amplification time is given by 



At w R-\/jIp/B r 



(12) 



and results to about 7500 yr for a star with R — 3R©, p = 
0.1 gem -3 , and an initial magnetic field strength of B r = 1 G. 
These time-scales are very short, but one needs to bear in mind 
that the stellar winds in this phase of the evolution will still remove 
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Figure 12. Dynamo-effect in terms of a^^ versus radius measured from the 
simulation with a perturbation at t = 0.003 = 9.5P ro t- The results from 
26 snapshots with 0.0035 < t < 0.004 of the simulation are superimposed 
in this graph. The top panel is for the northern hemisphere, the bottom panel 
is for the southern hemisphere. 



angular momentum and partly sustain the internal differential ro- 
tation. It is very likely that the whole process of amplification and 
onset of the instability take much longer in reality than suggested 
by our initial-value simulations. 

Instability occurs fairly early at an Alvfen angular velocity of 
about Qa = RB^ — 1000. This is much smaller than the condition 
for stability, Q 2 < derived by, e.g., Pitts & Tayler (1985). 
The difference is the much more complex structure of B^, in our 
case with much stronger currents than the configuration used for 
the analytical study, where B$ = s, B z — const. The fraction of 
the critical Alfven angular velocity to the actual rotation, Aa/O 
was even smaller in the stability analysis of the solar tachocline 
by Arlt, Sule & Riidiger (2007) where the toroidal magnetic fields 
were even more localised. 

Test runs on the linear stability (step II) indicate the a positive 
shear has a stronger stabilizing effect than negative shear in our par- 
ticular configuration. The test was performed on the toroidal field 
at t — 0.002 from step I with Bq — 300 and an arbitrary angular 
velocity fl(s) increasing by 15 per cent from s = to s = 1. The 
stability limit was at S > 2, while the stability limit for the same 
non-axisymmetric mode with dD./ds < was S = 1.2. This ex- 
plains the steep increase of the lines of marginal stability in Fig. [2] 
after t ~ 0.06. In the pre-main-sequence evolution, it may be re- 
sponsible for a long-term of stabilization as long as the star is gain- 
ing angular momentum on the surface by accretion. This situation 
is also stable against axisymmetric perturbations exciting the mag- 
netorotational instability. It will be very interesting to combine the 



simulations in this paper with the actual angular-momentum evo- 
lution obtained from combining the effects of accretion, magnetic 
star-disk coupling, and magnetized winds. 

During most of the pre-main-sequence period, the star is rotat- 
ing rather rapidly, with rotation periods of 1-2 days. Rapid rotation 
will suppress the Tayler instability in general. Once the accretion 
ceases and the disk turns into a more passive environment of the 
star, the stellar wind leads to a spin-down of the star. The increase 
of the rotation period can be quite steep, but only if strong magnetic 
fields are present (Stepien & Landstreet 2002). We are left with a 
contradiction here, since the spin-down requires poloidal fields, but 
the Tayler-instability providing these poloidal fields will set in only 
if the rotation period has already increased. 

Another question concerns the time it takes for the unstable 
magnetic field to become visible on the surface. In the 3D simula- 
tion, we find a maximum of Bmax = 168 at the surface, about 
18 rotation periods after the onset of the instability. Note that the 
initial poloidal magnetic field was internal; the surface field was 
zero, but it develops surface fields of the order of Bmax*' ~ 10 by 
diffusion. The bottom panel in Fig. [10] demonstrates that the mo- 
tions arising from the instability are mostly horizontal. The emer- 
gence of the flux at the stellar surface is therefore rather a diffusive 
process. In the numerical simulation, the dynamical and diffusive 
time-scales are not as widely separated as in reality. We may ar- 
gues that the emergence time is a diffusive one. Given a length- 
scale of say O.IRq for the rise and a microscopic diffusivity of 
1000 cm 2 s _1 , this time-scale results in about 9 Myr. The true rise 
time is certainly a combination of several effects and may be faster. 
The rise of a stable toroidal flux tube was studied by Mestel & 
Moss (2010). A tentative conversion of their time-scales into stel- 
lar values leads to 150 Myr and more which is more than an order 
of magnitude longer than in the unstable case (smaller structures). 

The maximum poloidal field at the surface converts to about 
2.7 kG using l[6}. The time-scale of 9 Myr is actually very inter- 
esting for Ap stars, it may take some time on the main sequence 
before the Ap phenomenon is observable (Hubrig, North & Mathys 
2000). Other observations favour a persistence of magnetic fields 
from early Herbig Ae/Be stars to main-sequence Ap stars (Wade 
et al. 2007). By contrast, Hubrig et al. (2009) find decreasing mag- 
netic fields in Herbig Ae/Be stars with age. While the observational 
picture will be gradually completing, simulations for a whole se- 
quence of magnetic Reynolds numbers will be need to tell whether 
the emergence time-scale is indeed a fraction of the diffusion time. 
Note that the above mentioned 2.7 kG are small compared to the 
internally possible toroidal fields of 10 5 -10 6 G. What appears to 
be a strong field in observations, can yet be a remnant of something 
much larger. 

Also the amount of complexity of the surface magnetic fields 
is compatible with surface fields from Zeeman Doppler imaging 
(e.g. Kochukhov et al. 2004). Note that differences in initial condi- 
tions may cause a variety of final surface field strengths as well as 
a variety of emergence times, very similar to the complex observa- 
tional picture that has been compiled up to now. 

We are left with two possible ways to magnetic Ap stars: the 
existence of fossil field configurations that are stable for very long 
times on the one hand, and the emergence of magnetic fields by 
an instability on the other hand, where the 'disrupted' configura- 
tion is the observed and long-term stable one. Our paper studies 
the latter option. Since the differential rotation vanishes due to the 
presence of the fields, there is no further build-up of toroidal fields 
and the 'disrupted' fields will be fairly stationary over evolution- 
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ary time-scales. A comparison of the topology at very late stages 
with the ones obtained in compressible, but non-rotating simula- 
tions by Braithwaite & Nordlund (2006) will be interesting, once 
long enough simulations are available. 

It is probable that the pre-main-sequence evolution of A stars 
is a mixture of various phenomena: the presence or absence of a 
substantial magnetic field in the collapsing cloud core, a possible 
short-lived dynamo during the convective phase (Arlt 2009), the 
angular-momentum evolution controlled by various torques, and 
the onset of the Tayler instability. It is still a challenge to explain 
why roughly 10 per cent of the possible scenarios lead to the Ap star 
phenomena, while the remaining fraction does not lead to strong 
fields at the stellar surfaces (fields may still be hidden in the star). 
However, the study of differences in the rotational evolution includ- 
ing differences in the amplitude of the differential rotation (deliv- 
ering both a wide range of toroidal fields and large difference in 
stability limits) are a promising field to find discriminating situa- 
tions between normal A stars and Ap stars. 

Dynamo action - Apart from looking at the the evolution of 
magnetic fields in the pre-main-sequence phase of Ap stars, we 
also studied the possibility of driving a dynamo from the instabil- 
ity in terms of representing it as a mean-field dynamo. We find 
enhanced values for the dynamo-Q and related effects describing 
the field generation in a mean-field context. An axisymmetric dy- 
namo driven by the Tayler instability cannot be excluded, although 
a proper energy source is missing in our setup. The prerequisite for 
the existence of an a-effect appears to be the presence of an inner 
boundary at which most of the a is concentrated, as is the kinetic 
helicity. The situation is thus interesting for Earth-like planetary in- 
teriors which have a solid inner core acting as an inner wall. To be 
interesting for a radiative zone, an inner boundary will be required 
such as a convective core, whose turbulent viscosity has a similar 
effect like a wall. This core should not be too far away from the 
Tayler-unstable zone, so a possible dynamo could be relevant for 
early B stars at best. The results do not support a dynamo from the 
Tayler instability in the solar tachocline though. 

An interesting domain is the Earth's dynamo though, since 
the outer, fluid spherical shell of the core is very similar to what 
was studied here. Because of the rather small differential rotation 
of the core, the onset of current-driven instabilities as well as he- 
licity distributions as shown in this Paper are possibly emerging in 
geodynamo simulations. 
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